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ABSTRACT 

The radial density profiles of stellar galaxy discs can be well approximated as an 
exponential. Compared to this canonical form, however, the profiles in the majority 
of disc galaxies show downward or upward breaks at large radii. Currently, there is no 
coherent explanation in a galaxy formation context of the radial profile per se, along 
with the two types of profile breaks. Using a set of controlled hydrodynamic simulations 
of disc galaxy formation, we find a correlation between the host halo’s initial angular 
momentum and the resulting radial profile of the stellar disc: galaxies that live in 
haloes with a low spin parameter A ^ 0.03 show an up-bending break in their disc 
density profiles, while galaxies in haloes of higher angular momentum show a down¬ 
bending break. We find that the case of pure exponential profiles (A « 0.035) coincides 
with the peak of the spin parameter distribution from cosmological simulations. Our 
simulations not only imply an explanation of the observed behaviours, but also suggest 
that the physical origin of this effect is related to the amount of radial redistribution 
of stellar mass, which is anti-correlated with A. 

Key words: hydrodynamics - methods: numerical - galaxies: spiral - galaxies: struc¬ 
ture. 


1 INTRODUCTION 


In his landmark work, Freeman (19701 found that most 


spiral galaxies share a uniform stellar surface brightness 
profile which is well fitted by an exponential, p* {R) oc 
exp (—Subsequent deeper imaging of a wider vari¬ 
ety of disc galaxies has shown some variation in the func¬ 
tional form of the radial profile. For example [Pohlen et al.| 
(20021 observed galaxies that exhibit two distinct exponen¬ 


tial proHles. In addition to an inner profile just as observed 


by Freeman (19701, they found a steeper outer exponen¬ 


tial profile such that the overall profile appeared to have a 


“break”. In barred SO-Sb galaxies Erwin et al. (20051 found 


profiles with an outer exponential slope which is shallower 


than the inner part. Pohlen & Trujillo (20061 compiled a 


sample of galaxies that included all three types of profiles. 
They called them pure, down-bending and up-bending ex- 
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ponentials. They found that only « 10% of their sample of 
late-type galaxies have pure exponential profiles that extend 
all the way out to the observational surface brightness limit. 
The measured abundances of down- and up-bending profiles 
are ~ 60% and ~ 30%, respectively. 


Using numerical simulations, Debattista et al. (20061 
explain down-bending disc breaks with stellar angular mo¬ 
mentum redistribution induced by the formation of a bar. 
A study by Roskar et al. (20081 that also employed numer¬ 
ical simulations predicted ‘U-shaped’ stellar age profiles for 
down-bending disc breaks and that the positions of the min¬ 
ima of the age profiles coincide with the break radius. This 


prediction was confirmed by observations (Yoachim et al. 

[Mo] ). 


There are many analytical models in the literature that 
study the physical origin of the exponential radial profile 


of stellar galaxy discs (e. g. Fall & Efstathiou 


1980 


Lin & 


Pringle||1987| [Dalcanton et al.||1997| |Mo et al.||1998[ 


den Bosch||2001 1. A common assumption of all these mod¬ 


els is that the distribution of specific angular momentum 
of baryons is conserved during the evolution of galaxies. 


Dutton (20091 found that low angular momentum material 
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needs to be removed in order to prevent centres of galaxies 
from becoming too dense. He suggested that stellar feed¬ 
back is a viable mechanism to redistribute angular momen¬ 
tum. Cosmological simulations of galaxy formation confirm 
this mechanism and find that this ejected low angular mo¬ 
mentum gas can be reaccreted with high angular momentum 


via the qalactic 

fountain effect (e. 

g. Brook et al. 12011 2012 

Marinacci et al. 

2011 

Ubler et al. 

2014 

I. However, currently 


there is no coherent explanation for the existence of pure, 
down- and up-bending radial profiles. 


In this Letter we present the first attempt at identifying 
physical quantities that determine the profile of stellar discs 
at large radii. We use a suite of smoothed particle hydro- 
dynamic (SPH) simulations of disc galaxies in an isolated 
set-up. We show that the type of disc profiles depends on 
the initial angular momentum of the galactic baryons. In 
Section we describe the simulation set-up followed by a 
presentation of the results in Section]^ and our conclusions 
in Section m 


2 SIMULATIONS 


In this Letter we present the results of a suite of simulations 
of disc galaxy formation from idealized and isolated, yet 
cosmologically motivated, initial conditions. We use a modi- 
fied version of the publicly available treeSPH c ode ChaNGa 
(Jetley et al. 2008 2010 Menon et al. 2014) ^ The simu¬ 


lations are evolved for 8 Gyr. This corresponds to 2 : ~ 1.5, 


when the last major merger era was coming to an end (Zent- 
ner fc Bullock|[2003 1. Examinations of cosmological simula¬ 
tions show that they behave similarly to isolated spheres af¬ 
ter their last major merger ( Zentner fc Bullock|2003 Obreja 
et ah, in preparation). In section 2.1 we present a detailed de¬ 


scription of our cosmologically motivated initial set-up while 
we will only briefly outline the implementation of hydrody¬ 
namics in section 1221 

Our initial set-up is motivated and indeed very simi¬ 
lar to that in [Roskar et al.| ( |2008P . We deliberately chose 
a simplified and controlled set-up in order to be able to 
link observed properties to their physical origin more eas¬ 
ily. This level of control comes at the expense of neglecting 
asymmetric influences from the cosmological context, such 
as merging of galaxies. However, stellar feedback helps to 
break the symmetry in the initial conditions. 


The gravitational softening is e = 227 pc. The SPH smooth¬ 
ing length is variable and set such that the kernel covers 50 
particles. 

The initial conditions were set up in four steps. First we 
created an equilibrium NFW DM halo ( Navarro et al.|1997 ) 
following the recipe from Kazantzidis et al. (20041 including 
an exponential cutoff outside R 20 Q ■ In the next step the mass 
of each DM particle was reduced by a factor of /b, the baryon 
fraction, and a gas particle was added at the same position, 
accounting for the mass difference between the old and the 
new DM particle. 

This gas sphere is then rotated by some random angle in 
order to prevent gas and DM particles from sharing identical 
positions. 

To set the gas velocities, we establish a cylindrical coor¬ 
dinate system (ur, Uc, Uz) such that the gas orbits about the 
z-axis. The velocities are set to obey the angular momentum 
profile for DM haloes as found W Bullock et al. (2001) in 
cosmological Abody simulation^ 


M (< j) 
M 200 


M- 


j/jn 


i/irnax -b/r - 1’ 


( 1 ) 


where M (< j) is the mass of all material that has less angu¬ 
lar momentum than j, fj. is the shape parameter and jmax is 
the maximum specific angular momentum in the halo, jmax 
depends on the value of /i and is proportional to the spin pa¬ 
rameter A. While A simply scales the gas particles’ angular 
momentum, /r sets the actual mass distribution of j/jmax. 
We use the definition of the spin parameter from |Bullock| 

eraLl([200Tl): 


A = 


J 


V2MVR 


( 2 ) 


^^=^^200 


Here J and M are the halo angular momentum and mass 
inside a sphere of radius R and V is the halo circular ve¬ 
locity at that radius. Radial and vertical velocities were set 
to Ur = Uz = 0. Tangential velocities are a function of the 
axisymmetric radius only (ut {R,(f),z) = Vt (R))- However, 
the DM halo does not rotate in our simulations. Finally, the 
gas temperatures were calculated such that the gas obeys 
hydrostatic equilibrium. 

Here we explore the effects of varying A at a fixed 
p = 1.3 on the radial profile of the resulting stellar disc. 
We explore the range 0.02 ^ A 0.][^ The simulation pa¬ 
rameters are summarized in Table [T] 


2.1 Initial conditions 

We set up isolated haloes with the following properties: 
M 20 Q = 1 O^^M 0 ,R 2 oo = 206kpc,/b = 0 . 1 ,c = 10 . M 200 
and R 200 are the virial mass and radius, ft is the fraction of 
baryons in the initial set-up and c is the halo concentration, 
/b is lower than the cosmological baryon fraction because our 
set-up does not account for high redshift outflows that are 
ejected from the halo. Dark matter (DM) particle masses are 

1.1 X 10® M 0 and initial gas particle masses are 1.2 x 10® M 0 . 


^ http://librarian.phys.Washington.edu/astro/index.php/ 
Research:ChaNGa 


2.2 Baryonic physics 


The ChaNGa code is derived from the treeSPH code gaso¬ 
line. It uses a modified version of SPH which employs 


a pressure averaged force calculation (Ritchie & Thomas 
|2001[ |Hopkin^|2013[ [Keller et al.||2014[). C haN Ga includes 
stochastic star formation (Stinson et al. 2006 c* = 0.1) 
based on a Kennicut-Schmidt law, radiative metal line cool¬ 
ing, metal diffusion and pre supernova stellar wind feedback 


^ These calculations were done assuming that the DM and the 
gas share the same angular momentum profile. 

® We did not explore lower values of A because the computational 
effort increases significantly as A decreases due to denser gaseous 
discs and an increased amount of star formation. 
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Table 1. Overview of all simulations and their properties. A is 
the initial spin parameter, Mgas and M* are the amount of gas 
and stars in the disc region {R < 30kpc, | 2 | < 3kpc) at 8 Gyr. 


A 

-Mgas 

(101° Mo) 

M* 

(101° Mo) 

0.02 

0.34 

3.14 

0.03 

0.56 

2.91 

0.035 

0.67 

2.64 

0.04 

0.95 

2.53 

0.045 

1.13 

2.42 

0.05 

1.24 

2.36 

0.055 

1.31 

2.30 

0.06 

1.34 

2.25 

0.1 

1.44 

1.96 


(early stellar feedback; Stinson et al. 20131. The feedback 


follows pDalla Vecchia fc Schaye (20121 in which the energy 
output from supernova explosions of a stellar population is 
released at one time altogether. The energy released per 
supernova is Ssn = 1.5 x 10®^ erg. Further details on the 
implemented physics in ChaNGa will be presented in up¬ 
coming papers | Keller et al. Stinson et al., in prepara¬ 

tion). First tests of this implementation produced realistic 
disc galaxies in cosmological simulations. 


3 RESULTS 


3.1 Disc profiles 

After evolving the simulations for 8 Gyr, we extracted the 
stellar surface density in axisymmetric radial bins E* (R). 
These stellar surface density profiles include all stars up to 
3 kpc above and below the plane and Inside a cylinder with a 
radius of 30 kpc. We fit these profiles with a superposition of 
a ‘broken exponential disc’ Ed surrounding an exponential 
bulge which is a very good parametrization of the data (see 
Fig.0: 

E* (R) = Ek exp f^ ^ -h Ed (R) (3) 


7?k 


where 


Ed (i?) = Eq X 




exp I 


Hb ) 
' RiJ 


exp 




ii R < Rh 
else 


(4) 

Here Ek and Eq are normalization factors for the bulge and 
disc component, respectively. Ri and Ro represent the Inner 
and outer disc scale-length and Rh the radius of the break. 
The probability distribution function (PDF) of the fit pa¬ 
rameters for the given surface density profiles was obtained 
using t he Monte Carlo Markov cha in (MCMC) algorithm 
EMCEE (Foreman-Mackey et al. 2013[ p| Fig .[^shows the stel¬ 
lar surface density profiles (black points) for all simulations 
at 8 Gyr overplotted by 100 models sampled from the ob¬ 
tained PDF (coloured lines). The grey data points show data 
for radial bins with nine or fewer star particles. These points 
were omitted in the fitting procedure. 


We did not use a standard minimization procedure because 
it gave unstable results. 


8 Gyr 



Figure 1. The radial stellar surface density profiles. Presented 
is the stellar surface density as a function of axisymmetric radius 
of the individual haloes at t = 8 Gyr. The errors are estimated 
as Poisson noise. The individual profiles are offset by 1 dex (i. e. 
one tick mark) each for clarity as we are not interested in the 
normalization. The grey data points correspond to radial bins 
with nine or less star particles which were ignored in the fitting 
procedure. The coloured lines show the model for 100 sets of 
parameters which were randomly chosen from the PDF obtained 
from the MCMC algorithm. For the up-bending disc profiles, the 
dashed lines are an extrapolation of the inner exponential part 
(Eq exp (—i?/Pi)) to make the break in the profile more easily 
visible. The figure qualitatively shows that there is a trend from 
up-bending disc profiles for galaxies in low spin haloes to down¬ 
bending profiles in high spin haloes. 


3.2 Disc breaks 

Fig.0 shows that the type of the disc profiles changes with 
A. In the lowest spin simulations (A ^ 0.03) the profile is up- 
bending, i. e. Ro > Ri- The models with high spin parame¬ 
ters (A ^ 0.045) clearly show a down-bending break. Fig. 
quantifies these trends in terms of the fit parameters. The 
middle panel shows the fitted inner to outer scale-lengths as 
a function of A. The inner scale-length Ri increases linearly 
with A. The outer scale-length Ro has a high value for low 
A and decreases with A until it approaches a constant lower 
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Figure 2. Systematic variation of the parameters for the bro¬ 
ken exponential disc profile (equation 13 ) as a function of A. 
The top panel shows the break radii. The estimates for the break 
radius for A = 0.035, 0.04 are plotted in grey since they are some¬ 
what ill-defined in these cases (see section [3^. The middle panel 
shows the estimated inner (black) and outer (red) scale radii. 
The bottom panel shows the ratio of the inner and outer scale- 
lengths. The horizontal dashed line indicates unity, that is the 
separation between the up- (below the line) and down-bending 
regime (above). The error bars indicate the range between the 
16-th and the 84-th percentile. The figure shows a clear trend 
from up- to down-bending profiles as A increases and a transition 
at A PS 0.035. 



R[kpc] 


Figure 3. Comparison of stellar surface density profiles relative 
to the current position of star particles (black) and their position 
at the time of formation (red). The top and bottom panels show 
the results for A = 0.03 and 0.06, respectively. The figure clearly 
shows that the amount of radial redistribution of stellar mass is 
much more pronounced for the case with small A. 


breaks (8-10 kpc) for other spin parameters. In such cases, 
the definition of _Rb and the distinction between Ri and Rc, 
are ill-defined. In the A = 0.04 case, the best fit for Rh is 
at the end of the disc. Therefore, the fitted break radius ex¬ 
ceeds the otherwise linear relation with A as can be seen in 
the middle panel of Fig.[^ In the A = 0.035 case, the best fit 
for Rh is at a minor wiggle in the profile. The result is that 
Rh does not fit the relation with A either but this time it 
falls short. This is not surprising given that the model pa¬ 
rameters are degenerate for pure exponential profiles where 
Rh and the distinction between Ri and Ro have no physical 
meaning. This phenomenon is reflected by the larger error 
bars. 


value for A > 0.04. The bottom panel of Fig. shows the 
ratios of the inner to outer scale-lengths Ri/Ro as a function 
of A. There is a transition from up-bending {Ri/Ro < 1) to 
down-bending {Ri/Ro > 1) disc breaks at A « 0.035. 

The top panel in Fig.j^shows the estimated break radii 
Rh as a function of A. Except for the A = 0.035, 0.04 cases, 
Rh shows a linear dependence on A. 

3.3 The pure exponential case 

In the region between the up- and down-bending regime 
(A = 0.035 — 0.04), the disc profiles are purely exponen¬ 
tial over a range that exceeds the typical position of disc 


3.4 Radial mass redistribution 

In this section we briefly lay out a possible mechanism that 
may be the cause for systematic changes of the disc pro¬ 
file with A. Fig. [^compares the final stellar surface density 
profile to the density profile of stars at birth, irrespective of 
their actual time of birth. Hereafter, we will refer to these 
as the final and formation profiles respectively. 

The two panels in Fig. show profiles for two different 
simulations, one with a spin parameter below the transi¬ 
tion region between up- and down-bending discs (A = 0.03, 
top panel) and one above (A = 0.06, bottom panel). In the 
A = 0.03 case we see that stars are formed with a very con¬ 
centrated profile that has a core in the centre. This core 
might be artificial due to centring issues. What is evident 
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without doubt is that practically all stars formed in the in¬ 
ner « 8 kpc of the disc, while the final profile extends all 
the way out to « 30 kpc. 

The figure shows that in the range 2 kpc ^ R ^ 5 kpc 
the formation profile significantly exceeds the final profile. 
Inside and outside that range the opposite is the case. This is 
because the integrated masses of the current and formation 
profiles are identical. It follows that a substantial amount of 
stellar redistribution has occurred from 2 kpc ^ R ^ 5 kpc 
outwards as well as inwards. 

For A = 0.06, we observe the same qualitative effect 
but there are two striking differences with important impli¬ 
cations. First, star formation extends to much larger radii 
(« 13 kpc). Secondly, despite the same qualitative trend the 
differences between the formation and final profiles are much 
smaller. The final profile closely follows the formation profile 
out to ~ 10 kpc and there is no central core in the formation 
profile. At large radii there is a steeply declining tail of stars 
that have reached large radii, leading to the down-bending 
profile. Therefore the amount of stellar redistribution that 
took place in this ‘high-spin’ case is much less than that 
in the ‘low-spin’ case. For clarity, we only presented these 
two cases. The intermediate cases show a monotonic trend 
between those shown in Fig. 


4 CONCLUSIONS 


Using numerical models of disc galaxy formation, we found 
a correlation between the initial spin of the host halo and 
the shape of the radial profile of the stellar disc. 

We find that galaxies with an initial spin parameter 
A < 0.035 form an up-bending disc profile while larger values 
of A yield down-bending discs. Pure exponential discs occur 
only right at A « 0.035 which coincides approximately with 
the median of the A distribution in cosmological simulations 
(A = 0.031; [Maccio et al.||2008| ). Thus, our model explains 
why only some late-type galaxies exhibit pure exponential 
disc profiles, while the majority of them show breaks in the 
outer disc profiles. As the transition between up- and down¬ 
bending disc profiles (A = 0.035) in our model coincides 
with the median of the A distribution, we expect roughly 
equal abundances of up- and down-bending profiles which 
are comparable but not equal to observational results (30 
% and 60 %; Pohlen fc Trujillo]|2006 K A possible cause for 
this discrepancy between our model and observations is that 
gas in cosmological simulations generally has larger angular 
momentum than the DM ( Stewart et al.|[2013 l. 

Nevertheless the model qualitatively reproduces all 
types of observed disc profiles. In a follow-up paper we will 
further explore this effect In cosmological zoom simulations. 
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